# 加载 fs 包
library(fs)


# 获取文件夹内的所有文件
file_list <- list.files(".\\HXcell_profiler", full.names = TRUE)

# 创建数据框
data_frame <- data.frame()

# 循环读取每个文件并将其加载到数据框中
for(i in 1:length(file_list)) {
  
  df <- read.csv(file_list[i], header = F)
  data_frame <- rbind(data_frame, df)
}
colnames(data_frame)<- gsub("V", "Feature", colnames(data_frame))

data_frame<-as.data.frame(t(data_frame))

colnames(data_frame)<-data_frame[1,]
data_frame<-data_frame[-1,]
colnames(data_frame)<-gsub("\\.", "-", colnames(data_frame))

library(Seurat)
library(patchwork)
library(reshape2)
library(RColorBrewer)
library(scales)
library(ggplot2)
library(ggpubr)
library(ggplotify)
library(pheatmap)
library(dplyr)
library(SeuratObject)
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor


meta<-read.csv("clusters_HX.csv",header = T,check.names = F)
meta$ID<- gsub("\\.", "-", meta$ID)
meta<-meta[meta$ID%in%colnames(data_frame),]

meta<-as.data.frame(t(meta))


colnames(meta)<-meta[1,]


meta <- meta[-1,]  # 这里我们假设你要替换的值为data_frame的第一列

meta<-as.data.frame(t(meta))

rownames(meta) <- gsub("\\.", "-", rownames(meta))

pbmc<-data_frame[,colnames(data_frame)%in%rownames(meta)] 




#pbmc <-pbmc [,1:30875]

pbmc <- CreateSeuratObject(counts = pbmc, project = "pbmc3k", min.cells = 3, min.features = 5)


# Normalizing the data
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

#Identification of highly variable features
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
# top10 <- head(VariableFeatures(pbmc), 10)
# plot variable features with and without labels
# plot1 <- VariableFeaturePlot(pbmc)
# plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
# plot1 + plot2

# Scale the data
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimPlot(pbmc, reduction = "pca")
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
ElbowPlot(pbmc)
pbmc <- FindNeighbors(pbmc, dims = 1:8)
pbmc <- FindClusters(pbmc, resolution = 0.17)
table(pbmc@meta.data$RNA_snn_res.0.17)

pbmc@meta.data$KNN<-meta$labels

pbmc@meta.data$KNN



pbmc <- RunUMAP(pbmc, dims = 1:8)

pA1 <- DimPlot(pbmc, reduction = "umap", label=T)+ NoLegend()


umap = pbmc@reductions$umap@cell.embeddings %>%  
  as.data.frame() %>% 
  cbind(cell_type = pbmc@meta.data$RNA_snn_res.0.17)

allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",  "#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",            "#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",            "#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
p <- ggplot(umap,aes(x= umap_1 , y = umap_2 ,color = cell_type)) +  geom_point(size = 1 , alpha =1 )  +  scale_color_manual(values = allcolour)

p2 <- p  +
  theme(panel.grid.major = element_blank(), #主网格线
        panel.grid.minor = element_blank(), #次网格线
        panel.border = element_blank(), #边框
        axis.title = element_blank(),  #轴标题
        axis.text = element_blank(), # 文本
        axis.ticks = element_blank(),
        panel.background = element_rect(fill = 'white'), #背景色
        plot.background=element_rect(fill="white"))
p2

p3 <- p2 +         
  theme(
    legend.title = element_blank(), #去掉legend.title 
    legend.key=element_rect(fill='white'), #
    legend.text = element_text(size=20), #设置legend标签的大小
    legend.key.size=unit(1,'cm') ) +  # 设置legend标签之间的大小
  guides(color = guide_legend(override.aes = list(size=5))) #设置legend中 点的大小 
p3

p4 <- p3 + 
  geom_segment(aes(x = min(umap$umap_1) , y = min(umap$umap_2) ,
                   xend = min(umap$umap_1) +3, yend = min(umap$umap_2) ),
               colour = "black", size=1,arrow = arrow(length = unit(0.3,"cm")))+ 
  geom_segment(aes(x = min(umap$umap_1)  , y = min(umap$umap_2)  ,
                   xend = min(umap$umap_1) , yend = min(umap$umap_2) + 3),
               colour = "black", size=1,arrow = arrow(length = unit(0.3,"cm"))) +
  annotate("text", x = min(umap$umap_1) +1.5, y = min(umap$umap_2) -1, label = "UMAP_1",
           color="black",size = 3, fontface="bold" ) + 
  annotate("text", x = min(umap$umap_1) -1, y = min(umap$umap_2) + 1.5, label = "UMAP_2",
           color="black",size = 3, fontface="bold" ,angle=90) 
p4


#pbmc <- RunTSNE(object = pbmc, dims.use = 1:8, do.fast = TRUE)
# V3

#DimPlot(pbmc,reduction = "tsne",label=T)



table(meta$labels)



pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()



library(magrittr)
library(monocle)
library(DDRTree)
library(ggsci)

#Extract data, phenotype data, and feature data from the SeuratObject
data <- as(as.matrix(pbmc@assays$RNA@layers$counts), 'sparseMatrix')
rownames(data)<-rownames(pbmc)
colnames(data)<-colnames(pbmc)

pd <- new('AnnotatedDataFrame', data = pbmc@meta.data)
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)

#构建S4对象，CellDataSet
HSMM <- newCellDataSet(data,
                       phenoData = pd,
                       featureData = fd,
                       lowerDetectionLimit = 1,
                       expressionFamily = negbinomial.size())



## 
HSMM <- estimateSizeFactors(HSMM)
HSMM <- estimateDispersions(HSMM)

HSMM <- detectGenes(HSMM, min_expr = 1 )
print(head(fData(HSMM)))

expressed_genes <- row.names(subset(fData(HSMM),
                                    num_cells_expressed >= 10))

head(pData(HSMM))
diff_test_res <- differentialGeneTest(HSMM[expressed_genes,],
                                      fullModelFormulaStr = "~ KNN")
ordering_genes <- row.names (subset(diff_test_res, qval= 1)) ## 不要也写0.1 ，而是要写0.01。

HSMM <- setOrderingFilter(HSMM, ordering_genes)
plot_ordering_genes(HSMM)

HSMM <- reduceDimension(HSMM, max_components = 3,
                        method = 'DDRTree') # DDRTree方式



HSMM <- orderCells(HSMM)

save(pbmc,meta,HSMM,
     file='immuneHX_data.Rdata')
rm(list=ls())
load(file='immuneHX_data.Rdata')

colour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",  
         "#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
         "#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
         "#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")

a1 <- plot_cell_trajectory(HSMM, color_by = "KNN") + scale_color_manual(values = colour)

a2 <- plot_cell_trajectory(HSMM, color_by = "seurat_clusters") + scale_color_manual(values = colour)

a3 <- plot_cell_trajectory(HSMM, color_by = "Pseudotime") 

p2 <- plot_complex_cell_trajectory(HSMM, x = 1, y = 2,
                                   color_by = "seurat_clusters")+
  scale_color_manual(values = colour) +
  theme(legend.title = element_blank()) 


CDS <- HSMM
CDS$KNN[CDS$KNN %in% c('1','2','6','0','7')] <- 'Other'
#CDS$KNN[CDS$KNN=='CS3' & CDS1$State %in% c(10:15)]<- '1A'
#CDS$KNN[CDS$KNN=='CS3' & CDS1$State %in% c(1:9,16:17)] <- '1B'
table(CDS$KNN)

color=c("#DC143C","#0000FF","#20B2AA","#F5F5F54D")

p1<-plot_cell_trajectory(CDS, color_by = "KNN") + scale_color_manual(values = color)


p2<-plot_complex_cell_trajectory(CDS, x = 1, y = 2,
                             color_by = "KNN")+
  scale_color_manual(values = color) +
  theme(legend.title = element_blank()) 

p1 / p2
#修改颜色透明度
library(grDevices)

color <- "#F5F5F5"
alpha <- 0.3

new_color <- adjustcolor(color, alpha)

# 输出修改后的颜色
print(new_color)
